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Abstract 


We propose the preconditioned global generalized product-type method 
based on the preconditioned global BiCG method to solve nonsymmet- 
ric saddle point problems with multiple right-hand sides. We apply an 
indefinite preconditioner to enhance the convergence rate of the method. 
We also present some theoretical analysis and discuss the convergence of 
the PGl-GPBiCG method. Some useful properties of the preconditioned 
matrix are established. Moreover, we present the bounds for the residual 
norm of the PGI-GPBiCG method according to the residual norm of the 
global GMRES method that guarantees convergence. Finally, some nu- 
merical examples are presented to show the efficiency of the new method 
in comparison with the preconditioned global BiCGSTAB method, and a 
comparison with another preconditioner is also provided. 
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1 Introduction 


Many problems in science and engineering applications such as mixed or 
mixed-hybrid finite element discretization of partial differential equations in 
computational fluid dynamics [13, 18, 31, 32, 41] and constrained optimiza- 
tion [29, 42, 43] eventuate in systems with multiple right-hand sides as 


jaro) be] = [rn e 


where A € R”*” is a symmetric positive definite matrix, B € R"*™ has a 
full column rank, X € R"**, Y € R™**, Fy € R"**%, and Fp € R™**. To 
convenience, we denote system (1) as follows: 


AX = B, (2) 


where ¥ = [¥M,..., 4] and B = [b,...,b“]. Under the assumptions 
mentioned above, the coefficient matrix A is nonsingular. 

Several efficient iterative methods have been proposed during the recent 
decades to solve the saddle point problems, such as the Uzawa method 
[11, 12, 45], the Hermitian and skew-Hermitian splitting (HSS) iteration 
methods [3, 5, 9, 35, 37], SOR-type schemes [7, 19, 20, 26, 40, 47] and 
Krylov subspace methods [1, 2, 10, 21, 33, 34], and so on. In order to im- 
prove the efficiency of standard iterative solvers, many preconditioners have 
been presented in the literature, for example, block diagonal preconditioners 
[28, 38], constraint preconditioners [6, 25], block triangular preconditioners 
[4, 14, 36, 44], parametrized block triangular preconditioners [24], and HSS 
preconditioners [5, 9, 22, 30]. 

In this paper, we study the preconditioned global BiCG (PGI-BiCG) 
method [23] and preconditioned global GPBiCG (PGI-GPBiCG) method [46] 
for solving the linear systems with multiple right-hand sides (1). We concen- 
trate on the use of the indefinite preconditioner 


ps esr 4 (3) 


This choice has been shown to be particularly effective on problems associated 
with constrained nonlinear programming [25, 27, 31]. In addition, in [33], 
the authors have shown that there is a tight connection between short term 
recurrence methods such as BiCG and the indefinite CG method used in [27]. 
More precisely, they are equivalent for a special choice of auxiliary vector, 
with which BiCG simplifies. In this paper, we discuss the convergence of the 
PGI-GPBiCG method and present two bounds for the residual norm of the 
method, which guarantee convergence. 

The outline of the paper is as follows. In Section 2, we mention some 
properties of indefinite preconditioner. In Section 3, we review the PGI] 
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BiCG method. In Section 4, we describe the PGl-GPBiCG method to solve 
the system (1). Section 5 is devoted to the convergence analysis of the global 
GMRES method. In Section 6, we study some theoretical properties of the 
convergence of the PGI-GPBiCG method. In Section 7, numerical experi- 
ments confirm the described theoretical results. In Section 8, we present our 
conclusions. 

Throughout this paper, we use the following notations. Inner product 
for two n x s matrices X and Y is defined as (X,Y) pr = tr(X7Y), where 
tr(Z) denotes the trace of the square matrix Z. The associated norm is the 
Frobenius norm denoted by || - || 7. We will use the notation (-,-)2 for the 
usual inner product in R”, and the related norm will be denoted by || - |l2. 
For a matrix V € R"*®, the block Krylov subspace K,(A,V) is defined by 
Kx(A, V) = span{V, AV, A?V,..., A*-1V}. Moreover, Z € K;(A,V) means 
that Z = )J=) & A’V, where €; € R, for i=0,...,j —1. Finally, 0, I,, and 
Ox will denote the zero, the identity, and zero matrices in R***, R***, and 
R!**, respectively. For brevity, we use the MATLAB-like notation [v; w] to 
represent the vector [v?w?]”. 


2 Properties of the indefinite preconditioner 


In this work, we use the global version of GPBiCG (Gl-GPBiCG) [46] for the 
solution of nonsymmetric saddle point problems with multiple right-hand 
sides (2). This method without a good preconditioner converges very slowly 
when applied to saddle point problems with multiple right-hand sides. In 
order to accelerate the convergence, we use the indefinite matrix P (defined 
in (3)) as a right preconditioner for the GL-GPBiCG algorithm applied to the 
problem (1): 


where 


n 


poe | I-III are) | 
— (BY By Be —+(BTB)"! ’ 
with IT = B(B7B)-1B", G = A(I—-Tl]) + Il, and S = 2(A— I)B(BTB)!. 
Once an approximate solution [X;; Ys] is determined, an approximate solu- 
tion to the unpreconditioned problem is recovered as [X%; Yk] = P~*[Xx; Yel. 
Choosing the vector [Xo; Yo] = [0; F2] as the starting approximate solution, 
the initial residual is given by 


me -am fal fe arn) 
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so that the second block component of Ro is identically zero. Problem (4) 
can thus be reformulated as determining an approximation [X;,; Y;] to the 
solution [X;Y] of the system 


Arp | = a (5) 
so that [Xz;Ye] = [Xo;Yo] + [Xn; Ye]. In addition, for every [U;0] € 


R™+")xs we have 


GJ] = Gl lolatul: 


In [33], the authors proved the following proposition for showing that the 
matrix G does have a full set of eigenvectors Z and for giving a bound for 
condition number of Z. 


Proposition 1. Let us assume that the matrix G = A(J — II) + II has 
nm —m nonunit eigenvalues \;,4 = 1,...,n —m. Let Z> be an orthonormal 
basis of span{II} and let the columns of Y; € R"*("—™ be eigenvectors of 
(I — I))AUW — ID) corresponding to all its nonzero eigenvalues. Then there 
exists an eigenvector matrix in the form Z = [Z,, Z| of A(J — II) + II such 
that 

I|All2 


< . i Som a, | 
K(Z)< (1+ lhl)" with ails fing Deal 


(7) 
where A = diag(A;) and y = ZF AY|(A—J)7?. 


The following proposition shows that the eigenvalues of G = A(J—II)+II 
belong to a real positive interval when A is a symmetric positive definite 
matrix. 


Proposition 2. Let 7 be an eigenvalue of G = A(J — II) + II. Then either 
7 =1or7T isa nonzero eigenvalue of (J — II)A(U — II). 


Proof. The proof is similar to that of Proposition 5 in [31]. 


3 The preconditioned global BiCG method 


In this section, we employ the PGI-BiCG method [23] to solve system (4). 
Let Xo € R”** be an initial guess with the residual Ryo = B—.A~'Xo and 
let Ro be an arbitrary n x s matrix. At step k, the residual R;, generated 
by this algorithm is such that, R, — Ro lies in the right matrix Krylov sub- 
space K;,(A,ARo) and R, is F-orthogonal to the left matrix Krylov subspace 
Kr (Ar Ro) = span{ Ro, AT Ro, sae AT ™ Ro}. 
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For solving (1) by preconditioned Global Biconjugate Gradient (PGI 
BiCG) algorithm by using the preconditioner P (defined in (3)) and equations 


(4) and (5), we choose [X; Yo] = (0; F)] as an initial guess. So, Ro = (75°), 


. ss = (1-1) RO ~ S 
As in [33], we set Ro = P~1Ro = (arp) aprR®)s Po = Ro, Po = Ro, and we 


I-IpR” ~ T-mpP 
ee and Py, = P*P. = ee, 
Therefore, the iterates R, and P; can be computed explicitly from R, and 
P,, and the auxiliary “tilde” recurrence can be omitted. Now, by using the 
relations (6) and ignoring from the last m rows of the matrices that are 
zero matrices, we can summarize the PGLBiCG algorithm for solving (2) as 
follows: 


obtain R, = P~1 Ry = ( 


Algorithm 1: The right PGI-BiCG algorithm for solving (1) 
Xo] _ [0 Xo] _ [0 
1. set [F°] = [2] ana [¥] = [1]: 
(1) 
2. Compute RS = F,—+(A-1)B(BT B)-1 Fy and set Ro = [Fs and Pp = 


Rw) 
0 


3. for k = 0,1,2,... until convergence 
<RO (-m RW >- 


4, = 
ue <@P® (7-0) POY sp 
5. Xyi = X_t+ a, P, Yet1 = Omxs 
1 1 1 2 
RO), = RO —a,GP™, RO), =0 
7 —_ <RUY), ,(I-M) RO > p 
<R® 1-1 RO sp 
@) _ pl) (1) (2) _ 
8. Pry = Rey + BePeo, Pray = 9 
9 end 7 - 7 7 7 
10. X41 = Xo + Xx41; _ Ye¢i = Yo+ Yesis 7 
Ll. Xpq1 = UF — W)Xp41 + 4 B(BTB)-1Y 41, Your = (BTB)-'BTXpa1 — 


(BT B)“*Y¥n44 


In practical implementation of Algorithm 1, we can factorize Bas B= QR 
and use the relation (B7B)~' = R7'R-?. We also note that the cost of 
solving with R is very low due to the particular structure and sparsity of the 
matrix. 

From Algorithm 1, the first block of residuals and first block of matrix 
directions can be expressed as follows: 


Rees = Re(G)RO, peas = Pr(Q)R, (8) 
where Rz(t), P.(t) € Px, Px is the set of polynomials p;,(t) of degree k with 


scalar coefficients satisfying p,(0) = 1. The polynomials R(t) and Px (t) are 
related together with the recurrence formulas as follows: 
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Rr+i(t) = Rr (t) — artPr (t), (9) 
Prti(t) = Reti(t) + BePr(t). (10) 


4 The preconditioned global GPBiCG method 


In the PGI-GPBiCG method [46], the matrix residual satisfies 
Ry = H(A 'P)Re(A'P)Ro, (11) 


where Ry is defined as above and Hz, is an accelerating scalar polynomial, 
which is computed as the following recurrence: 


Holt) =1, Golt) =, 
Hy (t) = He-1(t) — tGx_i(t), (12) 
G(t) = CeHe(t) + meGe_ift), k= 1,2,.... 


(1) 
If the PGL-GPBiCG algorithm is written for system (5) with Rp = [¥s i 


then the relation (11) can be written as 


1 
Ry = [He(G)Re(@)R”] 
0 
and for solving (1.2), the PGI-GPBiCG algorithm can be summarized as 
Algorithm 2. 


5 Convergence analysis of the global GMRES method 


In this section, we consider the block linear system GY = C’, where G € R"*” 
and Y,C’ € R”**, and we recall some convergence properties of the global 
GMRES method, which is needed in what follows. We consider the case 
where the matrix G is diagonalizable. Let G = 7DZ~—', where D is a diagonal 
matrix whose elements are the eigenvalues \1,..., An, and Z is the eigenvector 
matrix. In [8], the following upper bounds for the Frobenius norm of the kth 
residual Ry = C — GY, of the global GMRES method was given. Here 
Ro = C — GYp denotes the residual corresponding to the initial solution Yo. 


Theorem 1. [8] Let P;, be the set of polynomials of degree less or equal than 
k, and let K2(Z) = ||Z|l2||Z~+||2. Then we have the following results: 


(max  |p(A)|); (14) 


Rr < Z)||R 
|Relle < K2(Z)|| Rolle oes 


min 
pePx;p(0)=1 
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Algorithm 2: The right PGI-GPBiCG algorithm 
Xo _ {0 Xo _ {0 
L. Set Ea z B and ¥ = RE 
(a) 1 RY 
2. Compute Rj’ = Fi — +(A— 1) B(BB)~'Fy and set Ro = ; . 
3. Set Po = RS, Ty = W_y = P_1 =U_1 =Onxs, B-1 =0 
4. for k = 0,1,2,... until convergence 
5. Ph = RY + Br—1(Pr—1 — Up-1) 
6 — <R® (r-m) RY > p 
<GPx,(I-M) RO > 
7 Ye =Te_1 — RY — agWe_1 + onGPR 
8. Tr = RO — onGPr 
9 Cy = eat e¥k> P< CT Te > p< Vie Te> p<GT YE > B 
k ~~ <G@T,,GT,>F<YkiYk>F—<Yu,GTh>F<GTKYR>F 
10 — <GTy,GTp>r<Yr,Te>e-<Y¥e,GTp>r<GTp,Te>r 
- 7k <GT GT > #<YiVk> 8 —<Yb CTD F KOT VED F 
: Ty Th 
(if k = 0, then m = 0,¢, = 2o7,*aha= 
11. U; = Ch.GPr + x(Te—1 — RO + By—1Ug—1) 
12.0 2, = Gal 1k Zk—1 — AUR 7 
13.0 Xp41 = X~ + apPet Zr, Yeqi = Omxs 
14. RE, = Te — Yk — xT R= 
15. By = 28 <R@), ,-myRo? > 
: BT Ce <RO I -mROs 
16. Wry =GT, + BGP 
17. end 7 7 7 
18. Xp41 = Xo + Xk41; 7 Yeq1 = Yo+ Yer, 7 
19. Xpq1 = UF — W)Xp41 + +B(BTB) YG, Year = (BT B)-!BT Xp - 


(BT B)*Y¥n41 


where S'p(G) is the set of eigenvalues of matrix G. 


Theorem 2. [8] Let the initial residual Ro be decomposed as Ro = Z8, 


where 6 is an n x s matrix whose columns are denoted by 6),..., 6). 
Then . 
Z 
Rel < le _ 
ey (Vii yPVi+1) ey 
where 
r 4 1 Agere 
D= diag 1 ?,---, > 16PP} and Ven =[ i: : 
i=1 i=l 1 , ee dM 
The coefficients Bo, eng Ce) are the components of the vector 3, and e; is 


the first unit vector of R*t!. 
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From the fact that the eigenvalues of the matrix G = A(J — II) + II are 
positive, the min-max problem in the bound of Theorem 1 can be explicitly 
solved. Denoting by Amin and Amax the (positive) extreme eigenvalues of the 
matrix G, the following asymptotic bound holds: 


VE-1 
VK + 


where k = 3 stands for the ratio of the extremal (real) eigenvalues of the 


matrix G = A(I — II) +II (and so it is not its condition number!) [34]. 


I|Relle < &2(Z)|| Roll r( Ww (15) 


6 Convergence behavior of PGl-GPBiCG 


In this section, we propose an upper bound for the residual norm of the 
PGI-GPBiCG method. First, we derive an upper bound for the Frobe- 
nius norm of the kth residual RP GI-BiCG of the preconditioned global 
BiCG method. Using the steps 6 and 8 of Algorithm 1 and defining 
RO =(RP RM... RY), where RM, i=0,...,k, are the first block 
of the PGI-BiCG residual matrices, we can easily show that 


= = x 1 = 
GRY = ROP, — ——RO ET, 


Qk-1 


with EP = [0,,...,05,J,] € R°X** and Th = Ly Ao Up, where 


I, I, —Bols 
—I, Is : 
Ly = <1, ) Up = 7 i ry 
ee Is Is —Br-ols 
—I, I, Ts 
A, = diaglaoI,,...,@-11s], and a; and B;, i = 0,1,...,4, are scalars 


obtained from the PGIBiCG algorithm. In addition, T, is an invert- 
ible tridiagonal block matrix such that EPT Ey = az_1/,, where EP = 
[Lasery O30) eee, 

By assuming that the matrix RY, is of full rank and Gey = 
(GO RO a: as in [39], an upper bound for the norm of the 
first block of the residual of PGI-GPBiCG can be stated in the following 
theorem. 


Theorem 3. Suppose that GRY = ROT, - = RW ET with ETT Ey = 
1) 


apis, and V = [Iks Oxsxs](RO2,)*. If the matrix RW, is of full rank, 
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then for any polynomial p(t) of degree not exceeding k with p(0) = 1, we 
have : 7 
RPS Ile < Nall || PC@)RS? lle, (16) 


where Ny = In —(R® — RY FTV with FP = [I.,Is,..., 1s] € R°X**. 


Proof. The proof is similar to that of Theorem 1 in [39]. 


Now, we attempt to obtain an upper bound for the residual norm of 
the PGI-GPBiCG method. As in [39], we define the parameters ¢)41 = ¢;, 
7j+1 = ;, and matrices 

1 


Aya = HAO Re re wee, (eaeey = EGER PPS, j=0,...,k. 
J 


Using the relations (12), the iterates H; and G; can be computed by the 
recurrence formulas: 


Fim ge. = pores 
Hj41 = H; — GGG;, (17) 
Gr Tog Os 9 Ay 2a Re 
Cj+1 
By assuming that G #£0,7 =1,...,4 +1, we have 
Hyj41,Gj41 € IKa(G. a rome) Moreover, we obtain 


Hy = ere, (18) 


Under the assumption that all the generated coefficients ¢; are not zero and 


the grade yz of RW? o'-20G ith respect to G is not less than k, the recurrence 
formulas (17) determine the matrices 


By= Bi, Hy, ese5 Fy and Gy = (G1, Ga,...,Ga), 


whose matrix columns H; and G,;, 7 = 1,...,k, are linear independent. 


As in [39], by defining yj;41 = GnvaGa for j = 1,2,...,k, the relations 
(17) can be written as follows: 


Gb = ieae. Ae Cw, (19) 


where Aj = diag((J,,...;¢,ls|; BF = Og,...,0,,1,| € R°*”*, 
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I; —yls 
I, 
Ly = fx, and Uy, = : (20) 
a —YKLs 
I; 
Combining two equations in (19), we get 
GH, = Ay S, —G, Henk; (21) 


with Sk = Ly A, OU; and ETS, * Ey = Cpls, where ET = [1,,05,---, 0s] € 
R°xks, 
By assuming that the matrix Hj,, is of full rank and that H. ay 41 = 


(AE Apa) | Ee is similar to Theorem 3, we can present the following 
theorem for the PGI-GPBiCG method. 


Theorem 4. Let GH; = H,5; — C,H EP with ETS "Ey = CLs, and 


let V = [Tks pases) ga If the matrix H+, is of full rank, then for any 
polynomial p(t) of degree not exceeding k with p(0) = 1, we have 


|R@re-ersice|| . < || Myllr || p(G)RE POS | Ip, (22) 


where My = In — (Hg — Hast FT)V with FT = [I.,Ie,..-,To] €RO***. 


Proof. The proof is similar to that of Theorem 1 in [39]. 


Now, based on the above observations, we state the following theorem for 
bounding the residual norm of the PGI-GPBiCG method. 


Theorem 5. With the notation of Theorems 3 and 4 for Ro = and 


mo") 
0 
assuming that G = A(J — II) +I is diagonalizable, the right preconditioned 
global GPBiCG residual satisfies 


= i ' 7 Vvi—1 
[Bees (2(Z))?ILNell ell Melle|Rolle 75)", (23) 


where k = 4* stands for the ratio of the extremal (real) eigenvalues of 


matrix G = A(f— Il) +IL 
In addition, with the notation of Theorem 2 and RY) = BZ, we have 


- - Z K-11 
[REGHGes cals < K9(Z)||Nell el| Mall - Zl aE re ve 
vie (Vie41DVe+1)~*e1) 
(24) 


where D and Ve41 are defined in Theorem 2. 
So, the convergence is guaranteed. 
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Proof. Since the inequality (22) holds for any polynomial p(t) of degree not 
exceeding & with p(0) = 1, the use of the residual polynomial of the global 
GMRES (p@!-@“FES(G)) method in (22) implies that 


2 . - _GM Dears: 
RPO GPBICG| Z \| Mell | pe GueEe GR: )PGI-BicG Ile 


v 1 -—Bi K— 
< Ka(Z)|\Wall || ROPE-PLE Ip (LELYF, 


The second inequality follows from (15). Now, from the fact that (16) also 
holds for any polynomial p(t) of degree not exceeding k with p(0) = 1, we 
get 


1 _ A Ay - = 1 K— 
RQ PO SPF < Ko(Z)||Nell ell Melle || po OM BSG) RE Lp (YEGt)*. 


(25) 
This together with the inequality (15) implies that 


1 ss A vy 5 K— 
Rp Pe ePCF |e < (He2(Z))? Nell ell Mellel RG? ie(St)*. (26) 
Since RY = 0, so we have the inequality (23). 
Finally, due to Theorem 2, the inequality (25) yields the inequality given 
in (24). 


7 Numerical results 


In this section, we illustrate the efficiency of the PGlGPBiCG method for 
solving system (1). All the numerical experiments were performed in MAT- 
LAB R2017b with PC-Intel(R) Core(TM) i7, CPU 3.60 GHz, 16.00 GB of 
RAM. In all the examples, the starting guess was taken to be zero. The 
stopping criterion 

|Relle 


<10°° 
| Rolle 


was used. 


Example 1. In this example, a set of 6 problems was taken from the Univer- 
sity of Florida Sparse Matrix Collection[17] as the matrix A. These matrices 
with their generic properties are given in Table 1. Also, we consider € = 1, 
B = rand(n,m), Fy = rand(n,s), and Fy = rand(m,s), where the func- 
tion rand(l,k) creates an | x k random matrix with coefficients uniformly 
distributed in [0, 1]. 


We compare the Frobenius norm of the residuals (RES), the number of 
iterations (Iter), and the CPU time in seconds (CPU) required for conver- 
gence of the PGI-GPBiCG and PGI-BiCGSTAB methods for s = 5,10, 20. 
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Matrix Property 

order nnz spd 
1 | gr_30_30 | 900 7744 yes 
2 | nos3 960 15844 yes 
3 | 1138_bus | 11388 4054 yes 
4 | besstm11 1473 1473 yes 
5 | Dubcoval | 16129 253009 yes 
6 | bodyy4 17546 121550 yes 


Table 1: Test problems information 


The results obtained by these algorithms are presented in Table 2. As can 
be seen, the methods are similar in terms of the Frobenius norm of residu- 
als (|| Rx ||~), the number of iterations of the PGI-GPBiCG method is less 
than that of the PGI-BiCGSTAB method and the CPU time required for 
convergence in the PGI-GPBiCG method is similar or lower than that in the 
PGI-BiCGSTAB method except, for examples, nos3 with s = 10, 20. 


PGI-BiCGSTAB PGI-GPBiCG 
Matrix \ s 5 10 20 5 10 20 
Iter 47 49 49 38 38 38 
gr_30_30 | CPU 0.02 0.04 0.07 0.02 0.03 0.06 
RES _ 1.1082c-07 3.4408e-07 _6.3886e-07 | 2.6191e-07 4,2295e-07 _ 6.1987e-07 
Iter 192 183 193 167 180 172 
nos3 CPU 0.16 0.17 0.23 0.16 0.24 0.30 
RES —_ 2.7939e-07 7.3530e-07  8.0708e-07 | 3.9410e-07 7.6350e-07 — 1.0669e-06 
Iter 8364 4011 5393 2793 2677 2543 
1138_bus | CPU 9.61 6.02 10.72 4.23 4.84 6.33 
RES —_ 4.8730e-06 6.4503e-06 9.1015e-06 | 3.2011e-06 5.1722e-06 1.0123e-05 
Iter 335 381 441 216 335 250 
besstm11 CPU 0.98 1.20 1.72 0.59 1.05 1.11 
RES — 4.0741e-07 2.0770e-07 4.6741e-07 | 4.0252e-07 5.4683e-07  2.0180e-07 
Iter 105 99 99 85 85 86 
Dubcoval | CPU 30.48 35.21 40.87 24.62 31.08 37.33 
RES — 1.2814e-06 2.1262e-06 3.1051e-06 | 1.5025e-06 2.1308¢-06 _ 3.2268e-06 
Iter 182 194 170 165 166 165 
bodyy4 CPU 68.18 80.66 80.75 62.54 69.50 79.63 
RES _ 1.5969e-06 2.1782e-06 3.2651e-06 | 1.6013e-06 2.3639e-06 — 3.1067e-06 


Table 2: Numerical results for Example 1 with s = 5, 10,20 


Example 2. In this example, we consider the Stokes equation [16, 15] as 


—vAu+vp= f, in Q 
V:iu=g, inQ 
u=0, ond 
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where 2 = (0,1) x (0,1) C R?, 00 is the boundary of Q, v is the viscosity 
scalar, and u and p denote the velocity and the pressure, respectively. By 
discretizing (27), we obtain the system of linear equations as 


[-or ol [>= [-a]) 


in which 


_ TETI+TEeI 0) 2q7 x 2q7 _ 
a: 0 ees. ae 


I@F 


2q?xq? 
a es 


where T and F are tridiagonal matrices given by 


1 
i - -tridiag(-1,2,-1) ER™!, F = 5 tridiag(—1,1,0) €R™*, 
and ® denotes the Kronecker product. Also, h = —} is the discretization 


qt+1 
mesh size. We set n = 2q? and m = q?. Hence, the total number of variables 
is n +m = 3q?. We choose the right-hand side such that the exact solution 
of saddle point problem is a matrix of ones. For this example, we test three 
different v’s, that is, vy = 0.01,0.1,1 and q = 16, 32. 


In Table 3, we list the Frobenius norm of the residuals (RES), the number 
of iterations (Iter), and the CPU time in seconds (CPU) required for the 
convergence of the PGI-GPBiCG and PGI-BiCGSTAB methods with v = 
1,0.1,0.01, g = 16,32, and s = 5. As expected, the number of iterations of 
PGI-GPBiCG is less than that of PGI-BiCGSTAB. The CPU time obtained 
for PGI-GPBiCG is smaller than the one for PGLBiCGSTAB except for 
vy = 0.01 and q= 16. 


PGI-BiCGSTAB PGI-GPBiCG 
v\q 16 32 16 32 
Iter 38 74 23 47 
0.01 CPU 0.03 2.16 0.03 1.97 
RES  1.7694e-06 6.4937e-06 | 2.3360¢-06 6.6436e-06 
Iter 70 222 44 80 
0.1 CPU 0.04 3.10 0.03 2.20 
RES = 9.3235e-07 1.8710e-05 | 4.2220e-06 2.3113e-05 
Iter 83 828 37 82 
1 CPU 0.05 6.70 0.04 2.21 
RES — 5.5359e-05 2.7881e-04 | 3.5071e-05 2.3870e-04 


Table 3: Numerical results for Example 2 with v = 0.01,0.1,1 


In Figure 1, we display the convergence history of the PGI-GPBiCG and 
PGI-BiCGSTAB algorithms for Stokes problem with s = 5, q = 16,32, and 
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vy = 0.01,0.1,1, respectively. In these figures, the horizontal axis is the 
number of iterations (iters) and the vertical axis is the logarithm of the 
Frobenius norm of residuals (log;g || Rx ||7). The results show that whatever 
v becomes lower, the methods are smoother. Moreover, the PGl-GPBiCG 
method is more effective and smoother than the PGI-BiCGSTAB method. 


|——PGI-GPBICG 
PGI-BICGSTAB 


-——PGI-GPBICG 
PGI-BICGSTAB 


u i" " 
= 0 ae 4 
£ & u 
2-4 B..0 
2 4 
3 2 
4 3 
4 4 
oO 10 20 30 40 50 60 70 80 90 oO 100 200 300 400 500 600 700 800 900 
iters iters 
(a) v=1, q=16 (b) v=1, q=32 
3 4 
ay -——PGi-GPBICG ——PGi-GPBICG 
7 \ PGI-BICGSTAB 3 PGI-BICGSTAB 
1 eee 2 
0 1 N ‘ 
a4 eo 
eo :. e A 
Fy NN 2 \ 
3 2 
4 3 LL 
6 4 3 
4 5 
7 4 
oO 10 20 30 40 50 60 70 20 40 60 80 100 120 140 160 180 200 
iters iters 
(c) v=0.1, q=16 (d) v=0.1, q = 32 
3 3 
|——PGI-GPBICG |——PGI-GPBICG 
2 % PGI-BICGSTAB 2 ‘PGI-BICGSTAB 
1 * 1 
% . 
o o \ 
4 =4 \ 
Pa Pa ia 
£2 B2 
% 7 
3 3 
4 \ 4 \ 
5 5 
4 


5 10 15 20 Pa} 30 36 40 
iters 


(e) vy =0.01 q=16 


10 20 30 40 50 60 70 80 
iters 


(f) v =0.01, q = 32 


Figure 1: Convergence history of the PGI-GPBiCG algorithm for Stokes problem with 
s = 5 and different values of v and gq. 


Finally, we compare the numerical results obtained by the Gl-GPBiCG 
method with the preconditioner (3) and the preconditioner given in [1], 
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A B 
Pe.a.Q = Ee | 
with Q = I. For q = 16,32,64, v = 1, a = 0.05,0.1,0.5,1,10, and s = 5 
the numerical results are given in Table 4. As expected, the preconditioner 


discussed in this paper needs less CPU time than the preconditioner given in 
[1] except for g = 64 and a = 0.5, 1, 10. 


PGI-GPBiCG PGI-GPBiCG 
(pre. in [{1]) (pre. (3)) 
q\a 0.05 01 05 1 To 
Iter 196 88 25 22 35 37 
16 CPU | 0.66 0.35 0.11 0.10 0.14 0.04 
RES | 2.7200e-04 6.6880e-05 4.2144e-05 7.7757e-05 6.2674e-05 | 3.5071e-05 
Iter 276 109 38 35 53 82 
32 CPU | 17.10 6.87 2.50 2.35 3.53 2.21 
RES | 0.0013 6.1728e-04 1.5556e-4 —-.2.5773e-04 = 1.3334e-04 | 2.3870e-04 
Iter 318 162 56 56 91 201 
64 CPU | 315.67 162.91 60.20 60.15 94.10 128.57 
RES | 0.0060 0.0017 0.0021 0.0014 2.5978e-04 | 0.0015 


Table 4: Numerical results for Example 2 with a = 0.05, 0.1, 0.5, 1,10 


8 Conclusion 


In this paper, we applied the GI-GPBiCG method to solve the nonsymmetric 
saddle point problems with multiple right-hand sides. By using the indefinite 
preconditioner (3), we accelerated the convergence rate of the method. Also, 
we studied some theoretical properties of the PGl-GPBiCG method and pre- 
sented two bounds for the residual norm of the method, which guarantee the 
convergence. As expected, the experimental results showed that the number 
of iterations (Iter) of the PGI-GPBiCG method is less than that of the PGL 
BiCGSTAB method. The CPU time (CPU) of the PGlGPBiCG method 
is sometimes more than that of the PGI-BiCGSTAB method because of the 
existence of parameters ¢; and n;. In addition, the PGl-GPBiCG method is 
more effective and smoother than the PGI-BiCGSTAB method. 
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